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ABSTRACT 

We investigate the nonlinear properties of a system introduced by Burridge and 
Knopoff to model the dynamics of earthquakes. We find that a two-block system in a 
completely homogeneous configuration presents a complex behavior characterized by the 
presence of periodic, quasiperiodic and chaotic orbits. We have found routes to chaos 
via two types of intermittencies and period doubling bifurcations. The sensitivity of the 
evolution to different initial conditions is quantified by calculating the largest Liapunov 
exponent. The dynamics of the model is governed by nondifferentiable hows. 

PACS numbers: 05.45+b, 91.45.Dh, 91.60Ba. 
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Earthquakes are catastrophic events that happen when tectonic plates move with 
respect each other. The movement of the plates occurs in an intermittent way. A short 
period of slip follows a long period of rest, and in its turn is followed again by a short 
slipping period. This kind of dynamics is called stick slip motion. 

In 1967 Burridge and Knopoff [1] introduced two simple models that mimic the dy- 
namics of earthquakes. The systems consist of blocks connected by springs. The set is 
pulled with constant velocity on a surface with friction. In one of these models only the 
first block is coupled to the driving mechanism. It has been called "the train model" [2]. 
In the other model, all the blocks are connected to the driving element [1,3, 4]. We call 
it here the BK model. It is observed that in a completely homogeneous configuration of 
their component elements, both models present stick slip dynamics, and the distribution 
of the slipping events follows a power-law, in qualitative agreement with the distribution 
observed in real earthquakes, that is, the Gutenberg- Richter law[5]. However, there are 
fundamental differences between the train model and the BK model. In the first one the 
power-law distributions are limited only by the size of the system, that is, it is observed 
the existence of self-organized criticality[6]. In the BK model this does not occur. The 
power-laws have a limited extent [3,4] and they seem to result from finite size effects[7]. 

The studies on the Burridge-Knopoff models have been mostly concentrated to big 
chains, since the primary aim of these studies have been the comparison of the statistical 
distributions of the slipping events with the distribution of real earthquakes. In recent 
publications, the dynamics of a small chain of the BK model was investigated. Huang and 
Turcotte have shown that if an asymmetry is introduced in the friction forces, then a system 
of two blocks can present chaotic behavior [8, 9]. More specifically, they considered the case 
in which the friction force in one block is different from the friction force in the other block. 
They have compared their results with the motion of tectonic plates in California and in 
Japan and observed that they seem to have the same kind of chaotic dynamics observed 
in the two-block system of the BK model[9]. To the best of our knowledge, chaos has 
not been seen in the two-block system of the homogeneous BK model if no asymmetry is 
introduced in it. A one-block system, both in the BK and train models, does not present 
chaos. In this case only a periodic behavior is observed. 

The motivation for studying these mechanical models for earthquakes is in part to 
investigate one important open question in seismology, that is, are the temporal and spatial 
complexity of earthquakes generated by the nonlinearities of the equations that govern 
these phenomena? or are they due to geometrical and built-in heterogeneities? Another 
motivation results from the fact that these models are described by flows that are not 
infinitely differentiable, like most of the systems studied in chaos theory. Little is known 
about the class of systems we consider here, both from the physical and mathematical 
points of view. 
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The aim of this paper is to investigate the dynamics of a small system in the other 
model introduced by Burridge and Knopoff, i.e., the train model. In this case the two 
blocks would also represent two coupled tectonic faults. We find that without any kind of 
asymmetry, chaotic behavior is ubiquitous in a two-block system. We observe a rich bifur- 
cation diagram characterized by periodic, quasiperiodic and chaotic trajectories. We find 
entrances into chaos via intermittencies of types I and II and period doubling bifurcations. 
The phenomenon of crisis is also seen. To quantitatively characterize the divergence (or 
convergence) of nearby trajectories we calculate the most important (namely, the largest) 
Liapunov exponent of the system. 

The model we study is shown schematically in Fig. 1. It is the train model where 
the number of blocks is two. Each block of the system has mass m and the springs in the 
model are characterized by an elastic constant k. The first block is pulled with constant 
velocity v and the friction force F between the blocks and the surface is a function of the 
instantaneous velocity of the block with reference to a characteristic velocity v c . 

The equations of motion for the first and second blocks are given respectively by 

mX 1 = k(X 2 -2X 1 + vt)-F(X 1 /v c ), (la) 

mX 2 = k(X x - X 2 ) - F(X 2 /v c ), (16) 

where Xj denotes the displacement of the block measured with respect to the position 
where the sum of the elastic forces in it is zero. These equations are applicable only when 
the respective block is moving and the sum of the elastic forces in the block is larger than 
the maximum force of static friction. If this condition is not met we simply have Xj = 0. 
If we write the friction force as F(Xj/v c ) = F Q(Xj/v c ) where $(0) = 1 and introduce 
the variables r = u) p t, = k/m, Uj = kXj/F , Eqs. (1) can be written in the following 
dimensionless form[2] 

U 1 = U 2 -2U 1 + vt-$(U 1 /v c ), (26) 

U 2 = U 1 -U 2 -<S>{U 2 /v c ), (26) 

with v = v/V , v c = v c /V , and V Q = F j\Jkm. Dots now denote differentiation with 
respect to r. In a system of a single block the quantity F /u p is the maximum displacement 
of the pulling spring before the block starts to move; in the absence of dynamical friction 
2n /uj p and V Q are respectively a characteristic period of oscillation of the block and the 
maximum velocity it attains. We find that this system is completely described by two 
dimensionless parameters, v and v c . The system is four-dimensional, since its evolution 
is completely specified by giving the initial positions and velocities of the two blocks. We 
use here the velocity weakening friction force given by [3] 

*$!•>.) = rg*®-, (3) 
1 + \U\/u c 
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which is a simple nonlinear function. The friction force is the only nonlinear element in 
this model. 

A possible solution for the motion of the system is when both blocks move with 
constant velocity, equal to the pulling velocity v. The solutions for the positions of the 
blocks found from Eq. (2) are in this case 

U{ = -2/(1 + vjv c ) + ur, (4a) 

U% = -3/(1 + u/u c ) + ut, (46) 

where the superscript e denotes equilibrium position. The stability of this solution can by 
investigated by calculating the eigenvalues of the Jacobian matrix of a transformed system 
in which the equations are first order ODE's. This can be done by introducing two extra 
variables, V\ = U\ and V 2 = U 2 , which allow us to write the equations of motion as 

tfi = V 1 

V 1 = U 2 -2U 1 + vt-<S>(V 1 /v c ) 

(5) 

U 2 = V 2 

V 2 = U X -U 2 - ®{V 2 /v c ) 
The Jacobian matrix of this system when V\ = V 2 = v has the eigenvalues given by 



A ± JB± 

A, = (6) 

where i = 1, ...,4, A = v c l{y c + ^) 2 and B± = —6 ± 2v^5 + A 2 . For the region where 
we concentrate our attention (0 < v ^ 1 and < v~ x ^ 1) we find that B+ and B- are 
always smaller than zero. This gives two pairs of complex eigenvalues with 



|A-,| = |A 2 |= \/^-^ (7a) 



|A 3 | = |A 4 | = y^ (76) 

Since ^3,4) > 1, this implies that the equilibrium point is unstable. In fact, it is a saddle 
point with two stable directions and two unstable directions. Therefore, if this solution is 
perturbed, the motion settles into another attractor. 

For a random initial condition, we observe that the stick slip dynamics is the typical 
situation for the dynamical evolution of our system. Depending on the used parameters, we 
can find quasi-periodic, periodic or chaotic behavior. Examples for these types of motion 
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are shown in Fig. 2. In all the numerical simulations shown here we start the system with 
the blocks at rest. The initial position for each block is Uj = 0, i.e., with the sum of the 
elastic forces in the block equal to zero. Thus, no randomness is introduced, not even in 
the initial positions or velocities of the blocks. A transient time of r = 3000 is discareded 
in all simulations. Figs. 2(a) and 2(b) show the system evolution in phase-space for the 
first and second blocks, respectively, in a case of quasi-periodic motion where v~ x = 0.03. 
In Figs. 2(c) and 2(d) we show the dynamical evolution of the blocks when v~ x = 0.05, 
which gives a periodic motion. Note that for this case, only the first block sticks to the 
surface. Finally, Figs. 2(e) and 2(f) show the evolution when v~ x = 0.39, which gives a 
chaotic orbit. In all cases of Fig. 2 we have taken the pulling velocity v = 0.1. It is clear 
that the characterization of the motion as chaotic or nonchaotic can be done only with a 
more careful study, by investigating the Liapunov exponents. This has been done and will 
be discussed in the next paragraphs. 

Note that we have plotted in the x-axis the position of the block with respect to its 
(unstable) equilibrium position Uj, since the evolution of each block occurs around its 
respective equilibrium point. As seen in the figures, there is a discontinuity in Uj when 
the block sticks to the surface. So, this flow is not infinitely differentiable, like most of the 
flows considered in the literature of dynamical systems. 

By varying the degree of nonlinearity, i.e., varying v c , and fixing the pulling velocity 
we investigate how the motion changes its character from a periodic to a nonperiodic case. 
This is usually done by plotting bifurcation diagrams. For our case, a Poincare section 
[10] will lower the dimensionality of the system from four to three. Even in this situation, 
the visualization of the bifurcation diagrams is still complicated. We have verified that 
for our system the motion of the center-of-mass gives a good description of the dynamics 
of the model. Thus, we study bifurcation diagram by investigating the evolution of the 
center-of-mass in a determined Poincare section. However, we do not expect that this 
approach would always be good for a system with very large dimensionality. 

The coordinate of the center-of-mass with respect to the equilibrium position, deter- 
mined by Eq. 4, is W = (U x - U{ + U 2 - U%) /2. We take Poincare section ofWaXW = 0. 
In this way, we reduced the dynamics to a study of a one-dimensional variable. We show 
in Fig. 3 a bifurcation diagram for v = 0.1. On the x-axis we have v~ x and on the y-axis 
we plot W at W = 0. In the diagram we can see windows of periodic motion, and regions 
where the motion is nonperiodic. We find regions with period doubling bifurcation route 
to chaos, but also, as described later, other routes to chaos are also present in our system. 
We have studied bifurcation diagrams for other pulling velocities, and verified that their 
topologies can be different from the case we show here. These studies will be published 
elsewhere [11]. 

The nature of the nonperiodic motion can be analyzed by studying the sensitivity of 
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the dynamics to different nearby initial conditions. If this sensitivity is found, we have by 
definition a chaotic motion. Otherwise, the dynamics will be quasiperiodic. It is clear that 
in a periodic motion two nearby initial conditions will converge to the same orbit (if they 
are in the same basin of attraction) . In Fig. 4 we show the separation of two nearby orbits 
of the chaotic trajectory shown in Fig. 2(e). After the transient period dies out, we add 
a small perturbation of 10 to the positions and velocities of the blocks. At the instant of 
perturbation we have U\ — Uf = —1.72 and Ui = 0.09. We can see from the figure that 
at this point the two orbits practically coincide with each other. After a quite short time 
of integration (r = 60) presented in the figure, the two orbits have diverged considerably 
from each other. 

The quantity generally used to characterize the divergence (or convergence) of nearby 
trajectories is the Liapunov exponent. If the largest Liapunov exponent A m of the system 
is positive, then we have a chaotic motion. The calculation of A m can be done in the 
following way: After the transient dies out, we give small perturbations to the orbit and 
verify how much it separates after a small time interval from the original unperturbed 
orbit. The perturbations are given in the direction of the largest orbit separation. The 
logarithm of the average orbit separation along the trajectory, divided by the time interval 
gives A m . This standard method of calculating the largest Liapunov exponent of a system 
was introduced in [12]. 

In Fig. 5 we plot A m for the bifurcation diagram shown in Fig. 3. When the nonlin- 
earity is small, that is, v~ x is small, the motion is found to be periodic and quasiperiodic, 
as expected. For a flow this implies A m = 0, and the distinction of a motion with a large 
period and a quasi-periodic orbit can be done by calculating the second largest Liapunov 
exponent. However, here we concentrate our attention only to A m . The big period-three 
window loses stability at v~ x = 0.146 and after this the system enters into chaos. For 
v~ x ^ 0.146 the orbit stays for a long period of time close to the old period-three trajec- 
tory, and then makes chaotic incursions into other regions. That is, here we see what is 
called intermittency road to chaos. The intermittency is of type II, according to the clas- 
sification given by Pomeau and Manneville [13]. It seems difficult to calculate explicitly 
the eigenvalues of the map at the Poincare sections and in this way verify how they cross 
the unit circle, characterizing the intermittency. However, by studying the third Poincare 
return map, we see the typical structure of this type of intermittency. Details about this 
will be reported in a future publication [11]. After the first entrance into chaos, there exist 
several windows of periodic motion, giving A m = 0, as Fig. 5 shows. 

Other rich phenomena are also observed in the bifurcation diagrams. Period doubling 
cascades are born in tangent bifurcations, and present the well known Feigenbaum's ex- 
ponents [14]. A big periodic window is born in a tangent bifurcation at v~ x = 0.270. For 
v~ x ^ 0.270 there is an entrance into chaos via intermittency of type I [13]. Other type 
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I intermittencies are seen at v~ x = 0.498 and v~ x = 0.640. At v~ x rs 0.540 we see a 
crisis-induced intermittency [15]. That is, for v~ x }Z 0.540 the attractor suddenly widens. 
It spends long stretches of time in the region to which it was confined for v~ x ^ 0.540. At 
the end of these long stretches the orbit bursts out of the old region and bounces around 
chaotically in the new enlarged region made available to it by the crisis. Another crisis 
induced intermittency is seen at v~ x ps 0.683. For larger values of v~ x practically only 
chaotic behavior is seen. 

In conclusion, we have studied the dynamics of a two-block system in a one-dimensional| 
model for earthquake. We have found a rich dynamics with periodic, quasiperiodic and 
chaotic motion. The system presents several routes to chaos, as intermittencies of type 
I and II and period doubling bifurcations. A question we are currently investigating is 
whether or not the system becomes more chaotic as the number of block increases. We 
have seen that this model with large number of blocks presents self-organized criticality [2] , 
and it has been claimed that self-organized critical systems are not chaotic[16]. However, 
this contradicts the common result that a system becomes more chaotic as its dimension- 
ality increases. If chaos is a general solution for this model and if this system gives a 
reasonable description of the irregular dynamics observed in real earthquakes, this could 
mean that in practice earthquakes are predictable only on a short time basis, due to the 
sensitivity of the evolution to different initial conditions. 
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FIGURE CAPTIONS 



Fig. 1. System studied, which consists of a chain of two blocks connected by linear springs. 
The blocks are on a flat surface and the first one is pulled with constant velocity. 
Between the surface and the blocks there is velocity weakening friction force. 

Fig. 2. Orbits in phase-space for each block when v = 0.1. In (a) and (b) we have v~ x = 0.03, 
which results in a quasi-period orbit. For (c) and (d) we have v~ x = 0.05, which 
gives in a periodic trajectory. Finally, in (e) and (f) we have v~ x = 0.39, giving a 
chaotic motion. The vertical dotted lines just indicate the equilibrium position for the 
respective block. 

Fig. 3. Bifurcation diagrams of the velocity of center-of-mass W on the surface of section 

W = as a function of v~ x with v = 0.1. 
Fig. 4. Separation of two nearby orbits for the chaotic trajectory shown in Fig. 2(e) where 

v = 0.1 and u~ l = 0.39. A perturbation of 10 4 is given to the positions and velocities 

of the blocks when U\ — Uf = —1.72 and Ui = 0.09. The integration time for this 

figure is r = 60, after discarding the transient period. 
Fig. 5. The largest Liapunov exponent corresponding to the bifurcation diagram shown in Fig. 

3, where v = 0.1. The calculation of A m is done for an integration time r = 15000, 

with time steps of At = 0.01 and perturbations to the position and velocities of the 

blocks equal to 10~ 5 . 
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